In [1]:
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
%matplotlib inline
import pandas as pd
import argparse
import ellc
import gatspy
from gatspy import datasets, periodic
In [2]:
## V-band
In [3]:
V = pd.read_csv('V1061GettysburgV.csv', sep=",", header=0)
#X = pd.read_csv('V1061.csv', sep=",", header=0 )
HJDV = V["HJD"]
magV = V["Delmag"]
print (V)
In [4]:
# phase folded at
# http://vizier.cfa.harvard.edu/viz-bin/vizExec/Vgraph?J/ApJ/640/1018/./table3&P=2.35
In [5]:
plt.plot (HJDV[1:], magV[1:], 'o')
plt.title ('Magnitude vs. Julian Date')
ymin=float(min(magV[1:]))
ymax=float(max(magV[1:]))
plt.ylim(ymax, ymin)
plt.xlabel('Julian Date (d)')
plt.ylabel ('Magnitude')
Out[5]:
In [6]:
Period = 2.34665473 #days
PhaseV = (HJDV % Period)/ Period
plt.plot (PhaseV, magV, '-')
plt.title ('Magnitude vs. Phase')
ymin=float(min(magV))
ymax=float(max(magV))
plt.ylim(ymax, ymin)
plt.xlabel('Phase (d)')
plt.ylabel ('Magnitude')
Out[6]:
import numpy import ellc import matplotlib.pyplot as plt import astropy from astropy import units,constants from astropy.io import ascii
def afromP(m1, m2, P): return (((P2.) constants.G (m1 + m2) / (4*numpy.pi2.))**(1./3.)).decompose().to(units.AU)
Vobs = ascii.read('Torres2006_V1061Cyg_V.txt')
period = 2.34665473 * units.day
t_zero = 2452015.9050 #this is the spectroscopic only result (and it fits better)
m1 = 1.282 units.solMass m2 = 0.9315 units.solMass q = m2/m1 # = 0.7266 r1 = 1.615 units.solRad r2 = 0.974 units.solRad Teff1 = 6180. units.Kelvin Teff2 = 5300. units.Kelvin logg1 = 4.129 logg2 = 4.430 L1 = (10.0.533) * units.solLum L2 = (10.-0.173) units.solLum L3 = (10.**-0.153) units.solLum Mv1 = 3.456 Mv2 = 5.457 distance = 166.9*units.parsec
incl = 88.5 r_1 = 0.1625 r_2 = 0.0970
I1 = 0.737 I2 = 0.119 I3 = 0.144 sbratio = 0.472 # = J sigV = 0.0082
light_3 = I3
a = afromP(m1, m2, period)
pvalue = period.to(units.day).value avalue = a.to(units.solRad).value
phase = numpy.array([(tt % pvalue)/pvalue for tt in Vobs['HJD']])
shape_1 = 'sphere'#'roche'#'poly1p5' shape_2 = 'sphere'#'roche'#'poly1p5'
M_H = 0.0 #this is apparently unknown -- see Torres Section 5 : "No accurate estimate of the metallicity is available."
ldy_V = ellc.ldy.LimbGravityDarkeningCoeffs('V') a1_1, a2_1, a3_1, a4_1, y_1 = ldy_V(Teff1.to(units.Kelvin).value, logg1, M_H) a1_2, a2_2, a3_2, a4_2, y_2 = ldy_V(Teff2.to(units.Kelvin).value, logg2, M_H)
ld_1 = 'claret' ldc_1 = [a1_1, a2_1, a3_1, a4_1] ld_2 = 'claret' ldc_2 = [a1_2, a2_2, a3_2, a4_2]
sphase = numpy.argsort(phase) ellc_phase = phase[sphase] ellc_t = Vobs['HJD'][sphase] flux = ellc.lc(ellc_t, t_zero=t_zero, period=pvalue, #a=avalue, q=q, ld_1=ld_1, ldc_1=ldc_1, ld_2=ld_2, ldc_2=ldc_2, radius_1=r_1, radius_2=r_2, incl=incl, sbratio=sbratio, shape_1=shape_1, shape_2=shape_2, light_3=light_3)
ellc_mag = -2.5*numpy.log10(flux)
numpy.random.seed(seed = 1111) ellc_mag_obs = [numpy.random.normal(loc=x, scale=sigV) for x in ellc_mag]
print ("period = ",period) print ("semi-major axis = ",a) print ("inclination = ",incl) print ("mass1 = ",m1) print ("radius1 = ",r1) print ("radius1/a = ",r_1) print ("luminosity1 = ",L1) print ("Teff1 = ",Teff1) print ("log(g)1 = ",logg1) print ("mass2 = ",m2) print ("radius2 = ",r2) print ("radius2/a = ",r_2) print ("luminosity2 = ",L2) print ("Teff2 = ",Teff2) print ("log(g)2 = ",logg2) print ("surface brightness ratio = ",sbratio) print ("mass ratio = ",q) print ("light_3 = ",light_3)
Vuse = numpy.where(phase < 0.4)[0] Vsubtract = numpy.mean(Vobs['mag'][Vuse]) Vobs_plot = Vobs['mag'] - Vsubtract
plt.plot(ellc_phase, ellc_mag_obs,'o', label='ellc + noise', mfc='none') plt.plot(phase, Vobs_plot,'.',label='observations') plt.plot(ellc_phase, ellc_mag, label='ellc')
plt.ylim(max(numpy.append(Vobs_plot, ellc_mag)), min(numpy.append(Vobs_plot, ellc_mag))) plt.xlabel("phase") plt.ylabel("relative V magnitude") plt.legend()
plt.show()
In [7]:
#lc_dat = np.loadtxt("NightfallCurve.dat")
#phase = lc_dat[:,1]
#flux = 10**(-0.4*lc_dat[:,4]) # V-band
#M_H == 0
#Teff1 == 6180
#Teff2 == 5300
#logg1 == 4.129
#logg2 == 4.430
ldy_V = ellc.ldy.LimbGravityDarkeningCoeffs('V')
a1_1, a2_1, a3_1, a4_1, y = ldy_V(6180, 4.129, 0)
a1_2, a2_2, a3_2, a4_2, y = ldy_V(5300, 4.430, 0)
# Note that star 1 is the star labelled Secondary in nightfall
Vt_zero = 2452526.632
period = 2.34665473
Vsbratio = 0.472
Vr_1 = 0.1625
Vr_2 = 0.0970
Vincl = 88.5
light_3 = 0.144
a = 9.682
q = 0.7266
rotfac_1 = 1
rotfac_2 = 1
shape_1 = 'sphere'
shape_2 = 'sphere'
gdc_1 = -0.33
gdc_2 = -0.41
vsini_1 = 36
vsini_2 = 20
ld_1 = 'claret'
Vldc_1 = [a1_1, a2_1, a3_1, a4_1] #[0.50] # V-band coeffs from nightfall (using -Dv flag)
ld_2 = 'claret'
Vldc_2 = [a1_2, a2_2, a3_2, a4_2] #[0.60] # V-band coeffs from nightfall (using -Dv flag)
ecc = 0 # e>0 does not work in nightfall!
om = 90
f_c = np.sqrt(ecc)*np.cos(om*np.pi/180.)
f_s = np.sqrt(ecc)*np.sin(om*np.pi/180.)
Vt = PhaseV*period
lcV = ellc.lc(Vt,t_zero=Vt_zero, period=period, a=a, q=q,
ld_1=ld_1,ldc_1=Vldc_1,ld_2=ld_2,ldc_2=Vldc_2,
radius_1=Vr_1, radius_2=Vr_2,incl=Vincl,sbratio=Vsbratio,
rotfac_1=rotfac_1, rotfac_2=rotfac_2, gdc_1=gdc_1, gdc_2=gdc_2,
shape_1=shape_1,shape_2=shape_2, light_3=light_3, grid_1='default',grid_2='default')
magnV = -2.5*np.log10(lcV)
fontsize=9
fig=plt.figure(1,figsize=(8,4))
fig=plt.figure(1)
plt.subplot(211)
plt.xlim([-0.25,0.75])
plt.ylim([0.0,1.1])
plt.plot(PhaseV, lcV, linestyle='--')
plt.plot(PhaseV,magnV, linestyle=':')
plt.xlabel("Time [d]",fontsize=fontsize)
plt.ylabel("Flux",fontsize=fontsize)
plt.tick_params(axis='both', labelsize=fontsize)
plt.tight_layout()
plt.show()
In [8]:
noiseV = np.random.normal(0.00, 0.0082, 738)
magnV_noise = magnV + noiseV
#mean = 0.
#sigma = 0.0082
#tmp = np.random.normal(loc = mean, scale = sigma, size = 50)
#magn_noise_aaron = [np.random.normal(loc = x, scale = sigma) for x in magn]
#print(magn_noise_aaron)
#print(noise)
#print(np.array([1,2,3]) + np.array([4]))
In [9]:
plt.plot(PhaseV - 0.355, magV - 0.2, '-', label = "V1061 data")
plt.plot(PhaseV, magnV, '--', label = "ellc simulation")
plt.plot (PhaseV, magnV_noise, ':', label = "generated noise")
#ax.plot(Phase - 0.355, mag - 0.2, color='c', linestyle = '-', label="V1061 data")
#ax.plot(Phase, magn, color='indigo', linestyle = '--', label="ellc simulation")
#ax.plot(Phase, magn_noise, color= 'darkslategray', linestyle = ':', label = "generated noise")
plt.title('Change in Magnitude vs. Phase' )
plt.xlim([-0.50,1.50])
plt.ylim([0.6,-0.1])
#ymin=float(min(mag))
#ymax=float(max(mag))
#plt.ylim(ymax, ymin)
plt.legend(loc = 'best')
plt.ylabel('Change in Magnitude')
plt.xlabel('Phase (d)')
plt.show()
In [10]:
## R-band
# phase-folded at
# http://vizier.cfa.harvard.edu/viz-bin/vizExec/Vgraph?J/ApJ/640/1018/./table6&P=2.35
In [11]:
R = pd.read_csv('V1061GettysburgR.csv', sep=",", header=0)
#X = pd.read_csv('V1061.csv', sep=",", header=0 )
HJDR = R["HJD"]
magR = R["Delmag"]
print (R)
In [12]:
plt.plot (HJDR[1:], magR[1:], 'o')
plt.title ('Magnitude vs. Julian Date')
ymin=float(min(magR[1:]))
ymax=float(max(magR[1:]))
plt.ylim(ymax, ymin)
plt.xlabel('Julian Date (d)')
plt.ylabel ('Magnitude')
Out[12]:
In [13]:
Period = 2.34665473 #days
PhaseR = (HJDR % Period)/ Period
plt.plot (PhaseR, magR, '-')
plt.title ('Magnitude vs. Phase')
ymin=float(min(magR))
ymax=float(max(magR))
plt.ylim(ymax, ymin)
plt.xlabel('Phase (d)')
plt.ylabel ('Magnitude')
Out[13]:
In [14]:
#M_H == 0
#Teff1 == 6180
#Teff2 == 5300
#logg1 == 4.129
#logg2 == 4.430
ldy_R = ellc.ldy.LimbGravityDarkeningCoeffs('R')
a1_1, a2_1, a3_1, a4_1, y = ldy_R(6180, 4.129, 0)
a1_2, a2_2, a3_2, a4_2, y = ldy_R(5300, 4.430, 0)
# Note that star 1 is the star labelled Secondary in nightfall
Rt_zero = 2452526.633
period = 2.34665473
Rsbratio = 0.525
Rr_1 = 0.1635
Rr_2 = 0.0971
Rincl = 88.5
light_3 = 0.152
a = 9.681
q = 0.7266
rotfac_1 = 1
rotfac_2 = 1
shape_1 = 'sphere'
shape_2 = 'sphere'
gdc_1 = -0.33
gdc_2 = -0.41
vsini_1 = 36
vsini_2 = 20
ld_1 = 'claret'
Rldc_1 = [a1_1, a2_1, a3_1, a4_1] #[0.48]
ld_2 = 'claret'
Rldc_2 = [a1_2, a2_2, a3_2, a4_2] #[0.54]
ecc = 0 # e>0 does not work in nightfall!
om = 90
f_c = np.sqrt(ecc)*np.cos(om*np.pi/180.)
f_s = np.sqrt(ecc)*np.sin(om*np.pi/180.)
Rt = PhaseR*period
lcR = ellc.lc(Rt,t_zero=Rt_zero, period=period, a=a, q=q,
ld_1=ld_1,ldc_1=Rldc_1,ld_2=ld_2,ldc_2=Rldc_2,
radius_1=Rr_1, radius_2=Rr_2,incl=Rincl, sbratio=Rsbratio,
# rotfac_1=rotfac_1, rotfac_2=rotfac_2,
# gdc_1=gdc_1, gdc_2=gdc_2,
shape_1=shape_1,shape_2=shape_2, light_3=light_3)
# grid_1='default',grid_2='default')
magnR = -2.5*np.log10(lcR)
fontsize=9
fig=plt.figure(1,figsize=(8,4))
fig=plt.figure(1)
plt.subplot(211)
plt.xlim([-0.25,0.75])
plt.ylim([0.0,1.1])
plt.plot(PhaseR, lcR, linestyle='--')
plt.plot(PhaseR, magnR, linestyle=':')
plt.xlabel("Time [d]",fontsize=fontsize)
plt.ylabel("Flux",fontsize=fontsize)
plt.tick_params(axis='both', labelsize=fontsize)
plt.tight_layout()
plt.show()
In [15]:
noiseR = np.random.normal(-0.0081, 0.0081, 739)
magnR_noise = magnR + noiseR
In [16]:
plt.plot(PhaseR - 0.355, magR - 0.025, '-', label = "V1061 data")
plt.plot(PhaseR, magnR, '--', label = "ellc simulation")
plt.plot (PhaseR, magnR_noise + 0.01, ':', label = "generated noise")
#ax.plot(Phase - 0.355, mag - 0.025, color='c', linestyle = '-,' label="V1061 data")
#ax.plot(Phase, magn, color='indigo', linestyle = '--', label="ellc simulation")
plt.title('Magnitude vs. Phase' )
#plt.legend()
#plt.xlim([-0.50,1.50])
#plt.ylim([1.1,-0.1])
plt.xlim([-0.50,1.25])
plt.ylim([1.1,-0.1])
#ymin=float(min(mag))
#ymax=float(max(mag))
#plt.ylim(ymax, ymin)
plt.ylabel('Magnitude')
plt.xlabel('Phase (d)')
plt.show()
In [17]:
## I-band
# phase-folded at
# http://vizier.cfa.harvard.edu/viz-bin/vizExec/Vgraph?J/ApJ/640/1018/./table7&P=2.35
In [18]:
I = pd.read_csv('V1061GettysburgI.csv', sep=",", header=0)
#X = pd.read_csv('V1061.csv', sep=",", header=0 )
HJDI = I["HJD"]
magI = I["Delmag"]
print (I)
In [19]:
plt.plot (HJDI[1:], magI[1:], 'o')
plt.title ('Magnitude vs. Julian Date')
ymin=float(min(magI[1:]))
ymax=float(max(magI[1:]))
plt.ylim(ymax, ymin)
plt.xlabel('Julian Date (d)')
plt.ylabel ('Magnitude')
Out[19]:
In [20]:
Period = 2.34665473 #days
PhaseI = (HJDI % Period)/ Period
plt.plot (PhaseI, magI, '-')
plt.title ('Magnitude vs. Phase')
ymin=float(min(magI))
ymax=float(max(magI))
plt.ylim(ymax, ymin)
plt.xlabel('Phase (d)')
plt.ylabel ('Magnitude')
Out[20]:
In [21]:
#M_H == 0
#Teff1 == 6180
#Teff2 == 5300
#logg1 == 4.129
#logg2 == 4.430
ldy_I = ellc.ldy.LimbGravityDarkeningCoeffs('I')
a1_1, a2_1, a3_1, a4_1, y = ldy_I(6180, 4.129, 0)
a1_2, a2_2, a3_2, a4_2, y = ldy_I(5300, 4.430, 0)
# Note that star 1 is the star labelled Secondary in nightfall
It_zero = 2452526.633
period = 2.34665473
Isbratio = 0.590
Ir_1 = 0.1648
Ir_2 = 0.0994
Iincl = 88.3
light_3 = 0.156
a = 9.681
q = 0.7266
rotfac_1 = 1
rotfac_2 = 1
shape_1 = 'sphere'
shape_2 = 'sphere'
gdc_1 = -0.33
gdc_2 = -0.41
vsini_1 = 36
vsini_2 = 20
ld_1 = 'claret'
Ildc_1 = [a1_1, a2_1, a3_1, a4_1] #[0.35]
ld_2 = 'claret'
Ildc_2 = [a1_2, a2_2, a3_2, a4_2] #[0.42]
ecc = 0 # e>0 does not work in nightfall!
om = 90
f_c = np.sqrt(ecc)*np.cos(om*np.pi/180.)
f_s = np.sqrt(ecc)*np.sin(om*np.pi/180.)
It = PhaseI*period
lcI = ellc.lc(It,t_zero=It_zero, period=period, a=a, q=q,
ld_1=ld_1,ldc_1=Ildc_1,ld_2=ld_2,ldc_2=Ildc_2,
radius_1=Ir_1, radius_2=Ir_2,incl=Iincl,sbratio=Isbratio,
rotfac_1=rotfac_1, rotfac_2=rotfac_2, gdc_1=gdc_1, gdc_2=gdc_2,
shape_1=shape_1,shape_2=shape_2, light_3=light_3, grid_1='default',grid_2='default')
magnI = -2.5*np.log10(lcI)
fontsize=9
fig=plt.figure(1,figsize=(8,4))
fig=plt.figure(1)
plt.subplot(211)
plt.xlim([-0.25,0.75])
plt.ylim([0.0,1.1])
plt.plot(PhaseI, lcI, linestyle='--')
plt.plot(PhaseI,magnI,linestyle=':')
plt.xlabel("Time [d]",fontsize=fontsize)
plt.ylabel("Flux",fontsize=fontsize)
plt.tick_params(axis='both', labelsize=fontsize)
plt.tight_layout()
plt.show()
In [22]:
noiseI = np.random.normal(-0.0071, 0.0071, 740)
magnI_noise = magnI + noiseI
In [23]:
plt.plot(PhaseI - 0.355, magI + 0.16, '-', label = "V1061 data")
plt.plot(PhaseI, magnI, '--', label = "ellc simulation")
plt.plot (PhaseI, magnI_noise + 0.01, ':', label = "generated noise")
#ax.plot(Phase - 0.355, mag - 0.025, color='c', linestyle = '-,' label="V1061 data")
#ax.plot(Phase, magn, color='indigo', linestyle = '--', label="ellc simulation")
plt.title('Magnitude vs. Phase' )
plt.legend()
plt.xlim([-0.50,1.50])
plt.ylim([1.1,-0.1])
#ymin=float(min(mag))
#ymax=float(max(mag))
#plt.ylim(ymax, ymin)
plt.ylabel('Magnitude')
plt.xlabel('Phase (d)')
plt.show()
In [24]:
#M_H == 0
#Teff1 == 6180
#Teff2 == 5300
#logg1 == 4.129
#logg2 == 4.430
ldy_B = ellc.ldy.LimbGravityDarkeningCoeffs('B')
a1_1, a2_1, a3_1, a4_1, y = ldy_B(6180, 4.129, 0)
a1_2, a2_2, a3_2, a4_2, y = ldy_B(5300, 4.430, 0)
# Note that star 1 is the star labelled Secondary in nightfall
Bt_zero = 2452526.633
period = 2.34665473
Bsbratio = 0.398
Br_1 = 0.1635
Br_2 = 0.0959
Bincl = 89.1
light_3 = 0.137
a = 9.681
q = 0.7266
rotfac_1 = 1
rotfac_2 = 1
shape_1 = 'sphere'
shape_2 = 'sphere'
gdc_1 = -0.33
gdc_2 = -0.41
vsini_1 = 36
vsini_2 = 20
ld_1 = 'claret'
Bldc_1 = [a1_1, a2_1, a3_1, a4_1] #[0.74]
ld_2 = 'claret'
Bldc_2 = [a1_2, a2_2, a3_2, a4_2] #[0.83]
ecc = 0 # e>0 does not work in nightfall!
om = 90
f_c = np.sqrt(ecc)*np.cos(om*np.pi/180.)
f_s = np.sqrt(ecc)*np.sin(om*np.pi/180.)
Vt = PhaseV*period
lcB = ellc.lc(Vt,t_zero=Bt_zero, period=period, a=a, q=q,
ld_1=ld_1,ldc_1=Bldc_1,ld_2=ld_2,ldc_2=Bldc_2,
radius_1=Br_1, radius_2=Br_2,incl=Bincl,sbratio=Bsbratio,
rotfac_1=rotfac_1, rotfac_2=rotfac_2, gdc_1=gdc_1, gdc_2=gdc_2,
shape_1=shape_1,shape_2=shape_2, light_3=light_3, grid_1='default',grid_2='default')
magnB = -2.5*np.log10(lcB)
fontsize=9
fig=plt.figure(1,figsize=(8,4))
fig=plt.figure(1)
plt.subplot(211)
plt.xlim([-0.25,0.75])
plt.ylim([0.0,1.1])
plt.plot(PhaseV, lcB, linestyle='--')
plt.plot(PhaseV,magnB, linestyle=':')
plt.xlabel("Time [d]",fontsize=fontsize)
plt.ylabel("Flux",fontsize=fontsize)
plt.tick_params(axis='both', labelsize=fontsize)
plt.tight_layout()
plt.show()
In [25]:
noiseB = np.random.normal(-0.0071, 0.0071, 738)
magnB_noise = magnB + noiseB
In [26]:
plt.plot(PhaseV, magnB, '--', label = "ellc simulation")
plt.plot (PhaseV, magnB_noise + 0.01, ':', label = "generated noise")
#ax.plot(Phase - 0.355, mag - 0.025, color='c', linestyle = '-,' label="V1061 data")
#ax.plot(Phase, magn, color='indigo', linestyle = '--', label="ellc simulation")
plt.title('Magnitude vs. Phase' )
plt.legend()
plt.xlim([-0.50,1.50])
plt.ylim([1.1,-0.1])
#ymin=float(min(mag))
#ymax=float(max(mag))
#plt.ylim(ymax, ymin)
plt.ylabel('Magnitude')
plt.xlabel('Phase (d)')
plt.show()
In [27]:
## Condensation of above phase-folded light curves by band, theoretical and real data
In [32]:
plt.plot(PhaseV, magnV, label = "V-band")
plt.plot(PhaseR, magnR, label = "R-band")
plt.plot(PhaseI, magnI, label = "I-band")
plt.plot(PhaseV, magnB, label = "B-band")
plt.title('Simulated Magnitude vs. Phase in Multiple Bands' )
plt.legend()
plt.xlim([-0.25,1.50])
plt.ylim([0.7,-0.1])
plt.ylabel('Magnitude')
plt.xlabel('Phase (d)')
plt.show()
plt.plot(PhaseV - 0.355, magV - 0.2, '-', label = "V1061 data")
plt.plot(PhaseV, magnV, '--', label = "ellc simulation")
plt.plot (PhaseV, magnV_noise, ':', label = "generated noise")
plt.title('V-band: Magnitude vs. Phase' )
plt.legend()
plt.xlim([-0.50,1.50])
plt.ylim([1.1,-0.1])
plt.ylabel('Magnitude')
plt.xlabel('Phase (d)')
plt.show()
plt.plot(PhaseR - 0.355, magR - 0.025, '-', label = "V1061 data")
plt.plot(PhaseR, magnR, '--', label = "ellc simulation")
plt.plot (PhaseR, magnR_noise + 0.01, ':', label = "generated noise")
plt.title('R-band: Magnitude vs. Phase' )
plt.legend()
plt.xlim([-0.50,1.50])
plt.ylim([1.1,-0.1])
plt.ylabel('Magnitude')
plt.xlabel('Phase (d)')
plt.show()
plt.plot(PhaseI - 0.355, magI + 0.16, '-', label = "V1061 data")
plt.plot(PhaseI, magnI, '--', label = "ellc simulation")
plt.plot (PhaseI, magnI_noise + 0.01, ':', label = "generated noise")
plt.title('I-band: Magnitude vs. Phase' )
plt.legend()
plt.xlim([-0.50,1.50])
plt.ylim([1.1,-0.1])
plt.ylabel('Magnitude')
plt.xlabel('Phase (d)')
plt.show()
plt.plot(PhaseV, magnB, '--', label = "ellc simulation")
plt.plot (PhaseV, magnB_noise + 0.01, ':', label = "generated noise")
plt.title('B-band: Magnitude vs. Phase' )
plt.legend()
plt.xlim([-0.50,1.50])
plt.ylim([1.1,-0.1])
plt.ylabel('Magnitude')
plt.xlabel('Phase (d)')
plt.show()
In [29]:
## Lomb-Scargle Periodogram with gatspy
In [30]:
filts = ['V', 'R', 'I', 'B']
periods = np.linspace[0.2, 5.5, 1000]
model = periodic.NaiveMultiband(BaseModel=periodic.LombScargleFast)
model.fit(t, mags, dy, filts)
P = model.scores(periods)